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We compare numerical evolutions performed with the BSSNOK formulation and a conformal 
decomposition of a Z4-like formulation of General Relativity. The important difference between 
the two formulations is that the Z4 formulation has a propagating Hamiltonian constraint, whereas 
BSSNOK has a zero-speed characteristic variable in the constraint subsystem. In spherical symmetry 
we evolve both puncture and neutron star initial data. We demonstrate that the propagating nature 
of the Z4 constraints leads to results that compare favorably with BSSNOK evolutions, especially 
when matter is present in the spacetime. From the point of view of implementation the new system 
is a simple modification of BSSNOK. 



I. INTRODUCTION 

Numerical evolutions of the Einstein equations are 
complicated by the constraints of the system. At the 
continuum level the constraint subsystem of a given for- 
mulation of General Relativity (GR) must close; if the 
constraints are satisfied in one time-slice, they will re- 
main so. There are two possible numerical approaches to 
the constraints of the system. The first is to perform a 
so-called constrained evolution, in which the constraints 
are explicitly resolved at every point in time. There has 
been some progress in this direction in 3-D numerical rel- 
ativity since the work of Bonazzola et al. [IH (see also 
[ID, [3] for mathematical analysis and recent develop- 
ments), but the fact that the constraint equations are 
elliptic in character is a limitation from the point of view 
of computational cost. The second, more common ap- 
proach is to perform a free evolution. For free evolution 
the constraints are explicitly solved only in the initial 
time-slice. Throughout an evolution the constraints will 
be violated by numerical error, but this error should con- 
verge away with resolution. Computational resources are 
however always limited and in some situations a large 
constraint violation can be responsible for a numerical 
evolution to fail. Empirically spacetimes which contain 
matter are especially likely to suffer from large constraint 
violation. 

A major research topic has been the development of 
methods to minimise constraint violation at finite reso- 
lution. An obvious example is the PDEs and numeri- 
cal analysis of the constraint subsystem of a given for- 
mulation. Formulations with a constraint system whose 
characteristic speeds are non-zero are preferable, since 
one may hope that constraint violations propagate out 
of the numerical domain as opposed to sitting on the 
grid and growing. A related issue is the construction 
of constraint preserving boundary conditions, which pre- 
vent large constraint violations from propagating in from 
the outer boundary of the computational domain. Here 
we consider the effect of propagating constraints in com- 
bination with a damping scheme. 

The idea behind the constraint damping method is 



to identify places in the main system to subtract con- 
straints, such that the implied constraint subsystem picks 
up terms that force a reduction in the size of any viola- 
tion. In fact the constraint damping scheme of Gundlach 
et al. i 3:5] was used in the first successful numerical simu- 
lation of the binary black hole merger to dramatic effect. 
Essentially the same scheme is now used in the state of 
the art pseudo-spectral code of the Caltech-Cornell col- 
laboration,see for example (3|| and in the finite difference 
codes [HI, [131 • All of these successful uses of constraint 
damping have been for the generalized harmonic formu- 
lation of the Einstein equations, which implicitly relies on 
gauge conditions that are in the principal part harmonic. 

On the other hand, many groups have been using the 
BSSNOK formulation 0-Q (hereafter BSSN for brevity) 
with the 1+log and T-driver gauge conditions to perform 
simulations of both matter and vacuum space-times, fol- 
lowing the approach of [38|, |39| . Unfortunately, there has 
not been a systematic comparison of the constraint viola- 
tion between the different approaches, but every piece of 
evidence indicates that the simulations give comparable, 
compatible results. The differences between the general- 
ized harmonic and puncture approaches prevent the two 
methods from being straightforwardly interchanged. For 
example it is not possible to evolve the standard punc- 
ture initial data with the generalized harmonic approach. 
BSSN on the other hand lacks a constraint damping 
scheme. In this work we demonstrate that it is possible 
to evolve spacetimes of interest with a conformal decom- 
position of the Z4 [J, Q formulation, which we denote 
Z4c (conformal Z4) . Furthermore we show that the con- 
straint damping scheme for Z4 can be applied effectively 
and straightforwardly to the modified system. We find 
the new system particularly efficient in the evolution of 
non- vacuum spacetimes. From the implementation point 
of view, the new system is a simple adjustment of the 
BSSN equations. 

In section [TT] the Z4c system is described in some de- 
tail. Section IIIII contains a brief review of the equations 
of relativistic hydrodynamics. Our numerical results are 
contained in section|lVl Finally we conclude in sectionlVl 

Through the paper we use geometrical units c = G = 1 , 
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numerical results are reported with Mq = Mbh = 1- 



II. THE Z4C SYSTEM 



and 

The evolution equations are then 



(13) 



Formulation: The (constraint damped) Z4 formula- 
tion [H-Q takes the 4-dimensional Einstein equations and 
replaces them by 



Rab + y a Z b + V b Z a = Sn(T ab - ^g ab T) 
+ Ki[t a Z b + t b Z a - (1 + K 2 )g a bt c Z c }, 



(1) 



where Z a is a 4- vector of constraints and t a is a time-like 
vector field. Solutions of Eq. (P) are also valid solutions 
of the Einstein equations when the constraints Z a vanish. 
From the PDEs point of view, the most important part of 
the constraint addition is that of the partial derivatives. 
The constraint damping terms with coefficients Ki are 
discussed in more detail in the following section. We 
3 + 1 decompose the system and discard non-damping 
non-principal modifications to the ADM equations. This 
has the undesirable effect of breaking the 4-covariance 
of the Z4 formulation, but on the other hand we will see 
that discarding the lower order terms allows the evolution 
equations to be written very similarly to BSSN. The time- 
evolution equations are 

<9t7,j = Cp^j - 2aK i3 , (2) 

d t K i3 = —DiDjCt + a[R t3 - 2K lk K k + K tJ K + 2d (i Z j} ] 



Padm ) 



2 Si 



d t Q = a[-H + d k Z k - (2 + k 2 )ki6] + 
dtZi = ctMi + a9.j — an\Zi + fP Zij, 



(3) 
(4) 
(5) 



where Zi is just the spatial projection of Z a and 6 = 
—n a Z a . We stress that these equations of motion are 
not identical to those of Z4, although they differ only 
in non-principal terms. The Hamilonian and momentum 
constraints are 



H = R- l\,,l\'t + K 2 - 167rp A 



0. 



(6) 

M l = Dj(K iJ - f 3 K) - 8irS l = 0. (7) 



We adopt the standard notation for the metric variables 
and decompose the stress-energy tensor as 



Padm — n, a n b T 
Si = -^ a n b T ab , 

Sij = 7ial] b T ah . 



(8) 

(9) 
(10) 



To write the system as similarly as possible as the BSSN 
formulation we now conformally decompose the variables, 
or in other words make the change of variables 



Hi = 7 3 7u, 



k = Y j K i: 



29, 



X = 7 
A 



(11) 



n -l-HK lJ ~- llJ K), (12) 



dtX= -XHK + 2B)- A/? 4 ], 



(14) 



darj = -2aA %3 + p% jik + - (15) 

d t K = —D l D t a + aiAijlv + ^(K + 26) 2 ] 

+ 47ra[5 + p ADM ] + ^K ti + am (1 - k 2 )Q (16) 
the trace-free extrinsic curvature evolves with 
dtAij = xl-DiDja + a{R l3 - ZnS^f 
+ a[{K + 2Q)A ij -2A k i A kj ] 

+ p k A lhk + A tk (3 k J - ^A l0 P k tk (17) 

and finally we have 

d t f l = -2A ij aj + 2a[r jk A jk - 2A 11 ln(x),j 

- l^{k + 26),, - torjVSj] + y k (i] jk 

+ \f 3 P% + - f d% + ~f d 1 /^, (18) 

where we write ^ k T i j k = fV- The variable evolves 
according to Eq. Q with the appropriate substitutions 
Eqns. (|22«25l) . As with BSSN we write, in the d t in- 
equations 



Rii = R X ii + R 



(19) 



(20) 



~ -^DiXDjX- -^Jij&xDiX, 

Rij = --jl lm iij,lm + lk(i\^\J) +fd fc f(y) fc + 

f m (2ff (l f j)fem +fLf fe y). (21) 

The complete set of constraints are, in terms of the 
evolved variables, given by 

2Z i = 7 tf fJ'-7* fc 7y )fc) (22) 



e. 



H = R- A^Aij + -{K + 26) 2 - 16^p ADM , (23) 

W = djA^ + r 3k A 3k - \f j dj{K + 26) 

-|i«(log X )j, (24) 
D = ln(det 7) = 0, T = f j A i:j = (25) 

The change of variables introduces the algebraic con- 
straints D and T. The numerical implementation explic- 
itly imposes these constraints, so the continuum system 
we evolve is equivalent to Eqns. ([2][5]). 0ne may obtain 
BSSN from Z4c by taking 6 in Eqns. (fT4ll2l3|) . 
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Hyperbolicity: Statements about the hyperbolicity of 
Z4c obviously follow directly from those about Z4 since 
the two formulations share the same principal part. For 
completeness we summarize the relevant results here. A 
first order in time, second order in space system is called 
strongly hyperbolic if it has a complete set of charac- 
teristic variables with real characteristic speeds. This 
definition is equivalent to the existence of a strongly hy- 
perbolic first order reduction of the system [l(], El ■ The 
Z4c system is strongly hyperbolic when coupled to the 
puncture gauge conditions 



dta = (3 z a,i — a 2 [IlK, 



(26) 
(27) 



in which we always choose [Ll = 2/a and /xs = 1. To 
construct the characteristic variables we perform a 2 + 1 
decomposition in space and construct the fully second 
order characteristic variables |40| . Defining fas — 7^ us 
and 



d = ~(d t 
a 



the scalar sector has characteristic variables 



U s ±i 
U s±l ' = 



U s±±fi S 



d a± y/JIEa tS , 

do1qq±"fqq,s, 
2 2 

fis,S ~F 9^ Ps,Sl 

a or [is 

9oj3 s ± -^V^sA.s 



a 2 [i L 



do a 



(28) 



(29) 
(30) 



(31) 



afi s 



6(1 -As) 
a/is 



V3 
2 



Mi". 



3-4/is 



(32) 



with speeds ±{-JJx£, 1,1, 2a//*s/3) respectively, where 
As(3 - 4/2 s ) 



^ (3Mi-4As)(l-As)' 
The vector sector has characteristic variables 



(33) 



U 



U A 



±1 



d Pa ± VJlsPa,, 

do^fsA ± 7sA, 



/t s a 



(34) 

iflsl3 A , s ± do/3 A ), (35) 



with speeds ±(-y//ts, 1)- Finally the tensor sector has 
characteristic variables 



UaB±1 = dojAB ± lAB.si 



(36) 



with speeds ±1. Note that the choice fls — 1 renders 
some of the characteristic variables singular. In that spe- 
cial case one should analyse the system independently; it 
is again strongly hyperbolic, but we omit the character- 
istic variables. Strong hyperbolicity is enough to guaran- 
tee well-posedness of the Cauchy problem. For the initial 
boundary value problem, the stricter notion of symmetric 
hyperbolicity is desirable. We do not consider it here. 

Constraint subsystem: The principal part of the Z4c 
constraint subsystem is given by 



do® 
d Zi 
d Q H 



1 

2 

M, 



H + diZ\ 

2diM\ 
1 



doM, ~ -~diH + d^d,Z t - dt&Zj. 



(37) 
(38) 
(39) 

(40) 



Each of the variables satisfies a fully second order wave 
equation. The characteristic variables are simply 



d e ± e, s 

d H ± , 



doZi ± Zi tS 
d M, ± Mi. 



(41) 
(42) 



each with speeds ±1. In the numerical experiments that 
follow the crucial property is the propagating nature of 
the Z4c subsystem. The constraint damping scheme adds 
another level of complexity. However, if we linearize the 
system around flat-space we trivially recover the Z4 con- 
straint system; the analysis of @ holds immediately. All 
but constant frequency constraint violating modes are 
damped. 

Comparison with the BSSN constraint system: For a 
complete discussion of the BSSN formulation we sug- 
gest 043- Here we focus only on the BSSN constraint 
subsystem in order to compare with that of the Z4c sys- 
tem. Assuming that the algebraic constraints are satis- 
fied, the principal part of the constraint system is given 
by 



doZ t ~ M h doH ~ -2d t M\ 



(43) 



doM t ~ hdiH + &dj Zi + ^d i 8 j Zj . (44) 



The characteristic variables are 

c s±x = z SjS + \h± \m s 
c s0 = \h + z s , s . 



(45) 
(46) 



with speeds (±1,0) in the scalar sector and 

C A ±i = Z A ±M Al (47) 

with speeds ±1 in the vector sector. 

There are two important differences between the two 
constraint subsystems. The first is that BSSN does not 
accept a natural constraint damping scheme on every 
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constraint. The second is the 0-speed characteristic vari- 
able in the BSSN system. In free-evolution numerical ap- 
plications, some violation of the constraints is inevitable. 
Since the numerical scheme will inherit properties from 
the continuum system, one expects errors in the BSSN 
Hamiltonian constraint to sit on the numerical grid and, 
potentially, grow. On the other hand one should expect 
Z4c Hamiltonian constraint violations to propagate. In 
combination with suitable boundary conditions and the 
damping scheme we expect that Z4 Hamiltonian con- 
straint violations will be easier to control than those of 
BSSN. The conclusions obtained from the above linear 
PDEs analysis may not carry over to fully non-linear evo- 
lutions. Numerical evolutions are required to test their 
validity. 



III. RELATIVISTIC HYDRODYNAMICS 

We assume that the matter is described by a perfect 
fluid stress-energy tensor: 



T ab = phu a u b + pg ab 



(48) 



where p is the rest-mass density, e is the specific internal 
energy, h = 1 + e + pj p is the specific enthalpy, p is the 
pressure, and u a is the 4- velocity (u a u a — —1) of the 
fluid. The total energy density is given by e = p(l + e). 

The General Relativistic HydroDynamics equations for 
the perfect fluid matter {ideal GRHD) are the local con- 
servation law for the energy-momentum tensor, the con- 
servation law for the baryon number and the Equation of 
State (EoS) of the fluid: 



V a T Qb = 0, 
V a (pu a ) = 0, 
P(p,e)=p. 



Following 23;] we rewrite Eqns. 
flux-conservative form: 



(49) 
(50) 
(51) 

(|49H50P in first-order, 



by introducing the conservative variables: 
Q = Vrt A S k , t}, 

where 

D = Wp, 



S k = W 2 phv k , 
t = (W 2 ph 



P 



D 



(52) 

(53) 

(54) 
(55) 
(56) 



The simple physical interpretation of these variables is 
that they represent the rest-mass density (D), the mo- 
mentum density (Sk) and an internal energy (r = padm— 
D) as viewed by Eulerian observers. Above v % is the fluid 
velocity measured by the Eulerian observer: 



W 
W 



a 



1 



a \u 



(57) 



and W is the Lorentz factor between the fluid frame 
and the Eulerian observer, W = 1/ \/l — v 2 , with v 2 — 
7y« s u J . The fluxes in Eq. (|52|) are: 



i) 



D v 



S k 



v 

while the source terms are: 

T ab (d a g ak - T s ab g 5k ) , 
a(T M d a \na-T ab T ab )} 

= V=g{o, 

'1 



a 
a 



■pS, 



a 

•i 
fc) 



(58) 



pv 



(59) 



(60) 



-j3 l P 3 d k "iij - ad k a 



T oo (pipi Kij _ piQ.^ + T oi (2fi Kl] - d.a) + l ' ! K,,) 



Above g = detg ab = —a 7, 7 = det7y and Eq. ([60)1 
can be obtained using Eq. @ in a way to eliminate 
time derivatives. Note that both the fluxes and the 
source terms depend also on the primitives variables 
w = {p,p, e, u 1 }, and the source terms do not depend 
on derivatives of T ab . The system in Eq. (|5"2")) is strongly 
hyperbolic provided that the EoS is causal (the sound 
speed is less than the speed of light) (23[. 



IV. NUMERICAL RESULTS 

We consider the following set of tests in spherical sym- 
metry: 

Flat spacetime Evolution of constraint violating initial 
data. This test shows the basic mechanism of con- 
straint progation and damping at work; 

Puncture spacetime Evolution of puncture initial 
data. Z4c can be successfully used for the simu- 
lation of black-hole spacetimes with the puncture 
gauge; 

Stable star Evolution of a stable compact star. Con- 
straint propagation and damping are particularly 
useful in the simulation of non-vacuum spacetimes 
in which zero speed modes related to the Hamilto- 
nian constraint sit on the grid and grow during the 
evolution; 

Migrating star Evolution of an unstable star which mi- 
grates towards a stable one. The main features of 
Z4c are demonstrated in this nonlinear, very dy- 
namical scenario; 
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Collapsing star Evolution of an unstable star which 
collapses to a black hole. This test involves all of 
the difficulties of a simulation in numerical rela- 
tivity: matter, nonlinear dynamics, formation of a 
singularity and black hole evolution. 

The tests are also performed with BSSN for comparison. 
We stress that the differences between the results are 
convergent features, and with sufficiently high resolution 
they can be made arbitrarily small. In the presentation 
of the results we focus on the differences in the behavior 
of the Hamiltonian constraint. There is not a significant 
difference between the dynamics of the momentum and 
Zi constraints of the two formulations. In Z4c evolutions 
the constraint violation is typically of the same order 
as that of the Hamiltonian constraint. 

In each simulation we take for simplicity the constraint 
damping parameters «2 = and n\ = k = {0,0.02}, for 
a Z4c undamped/damped respectively. There is no par- 
ticular justification for the value k — 0.02. The value was 
found to be reasonable after the first numerical experi- 
ments. More detailed discussion follows. A systematic 
study of the k's is left for future work. 

The gauge choice is given by Eq. (|2"o) and Eq. (|2"7|) . one 
flavour of the popular puncture gauge [Tf|. In Eq. (|2"7)l 
we use rj = 2/Madm; choices more popular for matter 
evolutions, e.g. r) = 0.3/Madm HH, were tested in some 
cases, but without significant differences. 

The equations solved in these tests are obtained by a 
faithful spherical reduction of the 3D Cartesian equations 
presented in Sec. |TTJ Details are given in Appendix [A] 
Standard numerical methods based on finite-differences 
were used to solve them. They are described in Ap- 
pendix [Bj While in ID we may afford much higher reso- 
lutions than those of 3D simulations, we restrict to values 
reasonable for 3D mesh-refinement-parallel codes (see for 
instance Tab. [XJ) . All the figures refer to the highest res- 
olution simulations. 

Initial stellar models are built using the polytropic EoS 
(see Eq. (|B3p ) with adiabatic index T = 2 and polytropic 
constant K = 100 widely used in literature, e.g. Ref. [18j . 
The stable star model used in simulations in Sec. IIV CI 
has central density p c = 1.280 x 10~ 3 , gravitational 
mass M = 1.400 and circumferential radius R = 9.586 
(isotropic coordinate radius t\r = 8.126). The unstable 
model used in Sections IIVDI and IIV El has central den- 
sity p c = 7.993 x 10 -3 , gravitational mass M = 1.448 
and circumferential radius R — 5.838 (isotropic coordi- 
nate radius tr = 4.268). The EoS employed during the 
evolution is the ideal gas: 

P(p,e) = (T-l)pe. (61) 



A. Flat-space test 

We begin by evolving a large constraint violating per- 
turbation on flat-space in standard spherical coordinates. 
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FIG. 1: Hamiltonian constraint violation in the flat spacetime 
test for BSSN and Z4c. The main panel shows the violation 
in space at the initial and final time of the simulation. The 
lower inset focuses on large radii at the final time, where a 
pulse is propagating outwards (see text for discussion). The 
upper inset shows the 2-norm of the constraint in time. 



The violation is constructed by adding 
S X = exp(-r 2 /10) 



(62) 



to the x variable. According to section |TT] we should 
find that the Z4c Hamiltonian constraint will propagate 
away from the origin, whereas the BSSN Hamiltonian 
constraint should split into two components, one station- 
ary and one propagating. Fig. (p} demonstrates exactly 
the expected behavior. The main panel shows the Hamil- 
tonian violation in space at the final simulated time for 
the different formulations. In the BSSN evolution the vi- 
olation near the origin grows (cf. the dotted line, which 
is t = 0), while for the Z4c formulation it has propa- 
gated away. The propagation is clear from the lower in- 
set, which shows the violation at large radius; the initial 
violation propagates out almost completely in a form of 
a pulse in case of Z4c, and only partially in case of BSSN 
(the pulse is smaller). Note in addition that the damping 
mechanism is working effectively in Z4c with k = 0.02, 
for which the pulse become progressively smaller in space. 
Finally in the upper inset the 2-norms are presented. The 
global behavior of the constraint violation is better for 
Z4c as time advances. 

Note that it is possible to see over-damping effects in 
the Z4c simulations. In fact we experimented with dif- 
ferent damping values and for example with k = 2 we 
found a large constraint-violating "tail" left behind the 
outgoing pulse near the origin. 

The effect of constraint damping has been demon- 
strated mathematically around flat-space for non- 
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FIG. 2: Hamiltonian constraint violation in the puncture test 
for BSSN and Z4c. The 2-norm of the constraint is plotted 
in time. The inset shows the relative difference between the 
norms of the two simulations. 



FIG. 3: Oscillations of the central energy density in time in 
the star test for BSSN and Z4c. The inset shows the same for 
the whole simulation. 



constant in space violations with a plane wave ansatz. It 
is not entirely clear when that analysis will cease to deter- 
mine the effectiveness of the constraint damping scheme. 



B. Puncture 



Away from the puncture we find that there are small dif- 
ferences in the Hamiltonian constraint violation. All of 
the simulations in Fig. ([2]) were performed with Sommer- 
feld boundary conditions. Incoming constraint violation 
from the boundary can be seen at t = 200 in the inset of 
Fig. ©. 



We evolve a single stationary puncture with a precol- 
lapsed lapse and initially vanishing shift [35j . Non-trivial 
evolution occurs due to the gauge adjustment near the 
puncture. We find little difference in the evolution with 
the two formulations, as well as in violation of the Hamil- 
tonian constraint. 

Gauge dynamics dominate the evolution until roughly 
t ~ 25. After that the simulations settle at the stationary 
1 + log trumpet solution |34j [. Direct comparison with 
3D BSSN simulations of the BAM code [35J shows near 
identical evolution of the metric fields. 

The constraint violation results are summarized in 
Fig. @ which shows the 2-norm of the Hamiltonian con- 
straint during the evolution. The initial gauge-transient 
and the stationary phase are clearly recognizable. During 
the evolution the norm of Z4c is generally slightly better, 
see for example the inset in Fig. @ in which the 2-norm 
computed outside of the horizon is plotted. 

Looking at the spatial violation of the Hamiltonian 
constraint, there is large violation near the puncture, 
where the numerical solutions are almost identical. The 
violation is unsurprising since we are finite differencing 
across an irregular solution. Constraint damping appears 
to have no effect in a neighborhood of the puncture. 



C. Stable star 

We evolve a stable spherical star to t on d = 8000 
(around 40 ms). As is standard in such simulations, trun- 
cation error causes the star to oscillate at its proper ra- 
dial frequency. At late times oscillations are eventually 
damped by the numerical viscosity and by the interac- 
tion with the artificial atmosphere; a small linear drift of 
the mean value is usually also observed, 0, [l9|, [25[ . The 
numerical error is constraint violating. In this context we 
expect our linear PDEs analysis to be a good guide to the 
behavior of the non-linear system, since the unperturbed 
background solution is static and the numerical error is 
a small perturbation. 

Figure ([3]) shows the radial oscillations of the central 
rest-mass density in time. At the same resolution the 
amplitude of the oscillations in the Z4c evolutions are 
significantly reduced with respect to BSSN. Additionally 
at late times they are completely suppressed, as it is clear 
from the inset, and a drift, which improves in the con- 
strained damped Z4c, can also be seen. The frequencies 
of the radial mode and its overtones agree in both BSSN 
and Z4c with the correct values, see e.g. [2(| for both 
perturbative and 3D numerical results. 
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FIG. 4: Hamiltonian constraint violation in the stable star 
test for BSSN and Z4c. The zero speed mode of the BSSN 
constraint subsystem causes a slow growth; such growth does 
not occur in the Z4c evolution. The Z4c constraint damping 
scheme reduces the size of the violation. 



FIG. 5: Hamiltonian constraint violation in the stable star 
test. The 2-norm of the constraint up to r — 15 is plotted in 
time. The inset shows a zoom at early times. 



The dynamics of the Hamiltonian constraint in the 
BSSN and Z4c evolutions are similar at early times. Af- 
ter the interpolation from our spectral initial data grid to 
the evolution grid, finite difference derivatives are used 
to evaluate the Hamiltonian constraint. This produces 
a relatively large violation in the initial data at the sur- 
face of the star. The violation in space at the end of 
the simulation is plotted in Fig. The dotted line in 
the figure is the initial violation. In both systems we 
find that this violation propagates away from the surface 
of the star. However in the BSSN evolution, the prop- 
agation leaves behind a stationary violation that grows 
linearly in time. This accumulation affect in BSSN is 
demonstrated in Figs. (J3H5]). 

At the outer boundary a large violation of the con- 
straint is clearly visible in the BSSN evolution. It is 
caused by the Sommerfeld boundary conditions used with 
BSSN. This violation is non-convergent but it is usually 
mitigated by evolving with outer boundary further out, 
where the Sommerfeld conditions fare much better. The 
zero speed mode of BSSN plays the important role of 
keeping this violation at the boundary, thus minimizing 
their effect on the dynamics. On the other hand we find 
that for Z4c with Sommerfeld conditions the constraint 
violation from the outer boundary propagates inside the 
numerical domain and further perturbs the star, chang- 
ing for instance the value of the central density and excit- 
ing very large amplitude oscillations. The issue is com- 
pletely solved by the use of maximally dissipative bound- 
ary conditions (see Appendix [B]), which were used in all 
the matter simulations and in all of the presented figures. 



One possibility for not finding the same problem in our 
puncture evolutions is that in that case the constraint vi- 
olation from the boundary is typically much smaller than 
the violation at the puncture and so does not significantly 
alter the dynamics. Also in the black hole evolution we 
evolve with a much larger outer boundary. 

In Fig. (0 we plot the 2-norm of the Hamiltonian con- 
straint in time. The linear growth of the BSSN Hamil- 
tonian constraint is visible when the norm is computed 
up to r < 15 instead of the full grid. The norm on 
the full grid is in fact dominated by the violation caused 
by the Sommerfeld boundary conditions (see also again 
Fig. (U)). The reason is that the violation inside the 
star converges. At lower resolutions however the linear 
growth of the Hamiltonian constraint is visible in the 
norm on the full grid. Note that although Fig. ([5]) shows 
that the norm of the Hamiltonian constraint in the Z4c 
evolution drops below the initial value, the reason is the 
violation at the surface of the star in the initial data. 
Fig. (Q} shows that pointwise, at the end of the evolution 
the constraint violation has not dropped below the initial 
value except in a neighborhood of the surface of the star. 

Our experiments indicate that the propagation of the 
constraints is more helpful in preventing growth than 
constraint damping. 

D. Migration 

For the migration test we evolve initial data in hy- 
drostatic equilibrium that is unstable against first order 
perturbations. Truncation error causes the star to mi- 
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FIG. 6: Oscillations of the central energy density in time in 
the migration test for BSSN and Z4c. 



— BSSN - - Z4 k=0 Z4 k=0.02 




2000 4000 6000 8000 10000 
Time 

FIG. 7: Hamiltonian constraint violation in the migration 
test. The 2-norm of the constraint is plotted in time. 



grate to a model in the stable branch of the equilibrium 
configuration space. The energy difference between the 
two models is converted into kinetic energy and gener- 
ates large nonlinear pulsations. With an ideal EoS the 
transition involves the formation of a shock wave between 
an inner core and an outer (lower density) mantle which 
dissipates kinetic energy into thermal energy. Further- 
more a small amount of mass is expelled from the star 



when the shock reaches the surface. The dynamics are 
described in detail in [3, [l8| . This test is important be- 
cause it probes the behavior of the two formulations in 
a genuinely non-trivial situation, where the applicability 
of the linear analysis is unclear. 

Figure ([5]) shows the central density during the simula- 
tion: the strong nonlinear oscillations are clearly visible 
and each is accompanied by a bounce of the core of the 
star which feeds the shock. At early times the dynamics 
are comparable with very small differences. At late times 
(see the inset) however we find some differences between 
both BSSN and Z4c both with and without the damp- 
ing scheme. We believe the simulations at this point are 
not quantitatively reliable because of the lack of accuracy 
and convergence, so possible physical differences are not 
distinguishable from numerical errors. 

A second simulation carried with polytropic EoS (see 
Appendix [Bj Eq. (|B3p ). which excludes a priori the pres- 
ence of shock heating [4l|, completely solves this issue. 
Even at late times the solutions are comparable and the 
results obtained with BSSN and Z4c with k — are 
basically the same. Oscillations in this case are a lit- 
tle damped only by the numerical dissipation. We fur- 
thermore observe that the data obtained with Z4c and 
k = 0.02 are affected by some additional damping. 

The Hamiltonian violation behavior is qualitatively the 
same as in the stable star test. In Fig. ([7]) the 2-norm 
is plotted: the linear growth in the BSSN evolution is 
clearly visible (note that the boundary in this case is 
further out, so the Sommerfeld boundary condition con- 
tributes relatively little to the global violation) while the 
Z4c evolutions show a propagating and decreasing viola- 
tion. As opposed to the stable star case, at late times 
t > 2000 the 2-norm for Z4c damped is larger than Z4c 
undamped, but this is again a numerical artifact due to 
the reasons discussed above. The simulations performed 
with the polytropic EoS result in constraint violation sim- 
ilar to that in Fig. ([5]). 



E. Collapse 

Our final test is the evolution of the unstable star, 
which we perturb by a small (momentum constraint vio- 
lating) function 

Sv r = -0.005 sin(^), (63) 
R 

to trigger collapse to a black hole. The solution space 
of the two formulations is not the same when the con- 
straints are not satisfied. However we are able to follow 
the collapse with both BSSN and Z4c, and do not find 
significant differences in the dynamics of the two sys- 
tems. This can be seen for instance in the upper panel 
of Fig. (O in which we plot the central lapse in time. 
During the collapse the rest-mass is conserved to within 
a 1% error even in our lowest resolution simulation. 
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FIG. 8: Collapse dynamics: (upper panel) central value of 
the lapse and (bottom panel) irreducible mass of the black 
hole are plotted in in time in the collapse test for BSSN and 
Z4c. The inset in the bottom panel shows a zoom of the 
irreducible mass at late time. 



FIG. 9: Hamiltonian constraint violation in the collapse test. 
The 2-norm of the constraint is plotted in time. The inset 
shows the 2-norm computed only outside the apparent hori- 
zon; notice the log-scale of the y-axis. 



The Hamiltonian constraint violation of the formula- 
tions as the apparent horizon forms (t ~ 50) is compara- 
ble Fig. ([S]), but afterwards we find once again the linear 
growth in time of the BSSN Hamiltonian constraint. Fi- 
nally the lower panel of Fig. © shows that after the 
collapse the irreducible mass of the black hole: 



M irr = 



16tt 



(64) 



where A a h is the area evaluated from the apparent hori- 
zon, drifts from the correct value by roughly 1% in 
our BSSN evolution. At our lowest resolution (see Ap- 
pendix [Bj this drift is a serious (~ 25% by t = 300) 
problem, but it converges away at second order. On the 
other hand, the drift does not occur in any of the Z4c 
tests. Even at the lowest resolution the Z4c mass drift 
is below 0.1%. We mention that in Ref. [2l|, in which a 
different value of rj in the gauge condition is used, this 
effect is not seen. 



V. CONCLUSION 

Despite their success in the evolution of binary black 
hole and neutron star systems, the solutions of the free- 
evolution schemes widely used in numerical relativity 
simulations have some constraint violation. Motivated 
by the observation that BSSN evolutions of matter space- 
times typically have much larger Hamiltonian than (the 



almost negligible) momentum constraint violation, we 
considered the PDE properties of the BSSN constraint 
subsystem. The constraint subsystem has a zero speed 
characteristic variable involving the Hamiltonian con- 
straint, which we argue is the cause of the large violation. 
To demonstrate this with the aim of further improving 
the formulation for free-evolution schemes, we compare 
BSSN evolutions with those of a new conformal decom- 
position of the Z4 formulation. The latter formulation, 
referred as Z4c, has three desirable properties. Firstly 
the convenient choice of variables enables the evolution 
of puncture initial data. Secondly it inherits the propa- 
gating constraints of undecomposed system. Finally the 
constraint damping scheme described in Q can be at- 
tached to the new system. 

In a large set of tests in spherical symmetry, performed 
with faithful spherical reductions of the two Cartesian 
systems, we have demonstrated that the BSSN constraint 
violation can indeed be avoided by using a system with 
propagating constraints. We focus on the differences in 
the Hamiltonian constraint violation because other vari- 
ables do not show significant differences. We find that 
the propagation of the constraints is more helpful than 
the constraint damping scheme in controlling the viola- 
tion. We also find that when using the Z4c formula- 
tion boundary conditions are typically more important 
than when using BSSN. In our tests we find that maxi- 
mally dissipativc conditions arc always sufficient to avoid 
large constraint violation entering the grid from the outer 
boundary in Z4c evolutions. 
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Specifically we perform five tests, and focus on con- 
straint violation. We always find qualitatively the same 
behavior (similar orders of magnitude) in the momentum 
constraints of the different systems. Therefore we focus 
on the behavior of the Hamiltonian constraint. Firstly 
with a constraint violating perturbation on flat-space we 
find that the Z4c (BSSN) Hamiltonian constraint does 
(not) propagate, and that the Z4c constraint damping 
scheme works. Next we evolve a single puncture and find 
only small differences in the evolutions. When evolving 
either a stable or migrating star we find that at typical 
resolutions the Z4c Hamiltonian constraint violation is 
three or four orders of magnitude lower than that with 
BSSN. The migration test is especially important because 
it involves long-term nonlinear evolution. Our final test 
is the evolution of a star collapsing to a black hole. We 
find that the two formulations give very similar dynam- 
ics, which is not guaranteed because our initial data is 
constraint violating. We find that after the collapse the 
BSSN Hamiltonian constraint grows linearly in time, and 
that the irreducible mass of the black hole is not per- 
fectly conserved at low resolutions. These problems are 
not present in the Z4c tests. 

In the work presented here we have not tuned the con- 
straint damping parameters, so an obvious step is to in- 
vestigate the effect of these parameters in more detail. 
However the most important follow up is to investigate 
whether or not similar results can be reproduced in 3D, 
where there may be additional issues affecting numerical 
stability. If 3D astrophysical simulations with Z4c are 
possible it will be interesting to see whether or not the 
improved Hamiltonian constraint behavior will effect the 
physics, in particular the extracted gravitational waves 
from binary systems or stellar collapse. 



in part by DFG grant SFB/Transregio 7 "Gravitational 
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Appendix A: Spherical reduction of the equations 

In spherical symmetry the line-element is: 

ds 2 = - (a 2 - I3 r p r ) At 2 + 2/3 r drdt + 7rr dr 2 +7 T r 2 dn 2 , 

.( A1 ) 

where f3 r = (3 r /"/ rr - The determinant of the 3-metric can 
be written as: 
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: 6>7 S 



(A2) 



with the definition 7 s = 7 rr 7^. For the spherical reduc- 
tion of the metric equations we follow 1_7J. An auxil- 
iary flat derivative whose connection vanishes in Carte- 
sian coordinates is defined. The standard formulation of 
the equations of motion is given in Cartesian coordinates. 
Thus the spherical reduction is made by replacing partial 
derivatives in constraint addition terms and non-tensorial 
variables with the auxiliary derivative, and simply rotat- 
ing from Cartesian to spherical coordinates assuming two 
angular killing vectors. The resulting expressions are too 
unwieldy to be presented here. 

The spherical reduction of the hydrodynamics equa- 
tions of motion (|52p is performed defining conservative 
variables as: 



q = V¥{D,S r ,T} 



(A3) 



The equations result in a flux-balance form with the 
fluxes given by: 
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Note that, as in the field equations, all the source terms are regular, due to the regularity conditions at the origin, 
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TABLE I: Settings used in the simulations. Note that dimen- 
sionless units c = G = Mq — Mbh = 1 are employed. 

Test Resolution, N Grid Atmosphere Time 

{low,med,hig} r out /thr, /lev tend, c c a 

Flat {1000,2000,4000} 100 — , — 70,0.5 

Puncture {2000,4000,8000} 200 , 400,0.5 

Stable Star {100,200,400} 20 10~ 8 , 10" 12 8000,0.5 

Migration {300,600,1200} 40 10" 6 , 10" 9 10000,0.25 

Collapse {400,800,1600} 20 10~ 7 , 10~ 12 300,0.5 



since 1/r terms appear always together with a vector 
quantity. Note also that the system is slightly different 
from a previous formulation, see e.g. Ref. |26|,[27j], and in 
our case the characteristic speeds are the same as the ID 
Cartesian one given in [23| , once the formal replacements 
747 s , v x ^ v r and r ) xx j rr are applied. 

Appendix B: Implementation 

Matter and metric helds are integrated forward in time 
together using the Method of Lines (MoL) based on 
Runge-Kutta (RK) inte gra tors. Both TVD 3rd-order and 
4th-order RK, see e.g. [28|], were used for the computa- 
tions subjected to the Courant condition for the timestep 
At = c c flAr. 

The RHS of the equations are discretized on a stag- 
gered grid in r £ (0, r out ] with uniform spacing Ar = 
fout/N . Spatial derivatives of the geometric helds are dis- 
cretized using centered 4th-order hnite differences, and 
the 4th-order Kreiss-Oliger dissipation operator is always 
used in simulations, with a — 0.007. The analysis of [33[ 
demonstrates that the standard discretization of the Z4 
formulation in terms of the ADM variables, coupled to 
the Bona-Masso gauge and vanishing shift, will result in 
an unstable scheme if artificial dissipation is not used. 
However It is not immediately clear whether or not that 
analysis holds in the conformally decomposed Z4c case. 
All the relevant parameters used are listed in Tab. Q] 

We use two different boundary conditions on the metric 
variables; Sommerfeld and maximally dissipative. With 
BSSN we always use Sommerfeld conditions like those 
typically used in 3D implementations. At the boundary 
each variable is modeled by flat space plus a perturbation. 
One may then derive the boundary conditions 

dtf = -vf,r--(f-fo), (Bl) 
r 

where /o is the background solution of the field. For Z4 
we also use maximally dissipative boundary conditions. 
To construct the conditions we linearize the spherical 
equations of motion around flat-space and analyze the 
system in fully second order form. We construct the in- 
coming characteristic variables and use conditions 

U- = 9. (B2) 



The given data g for the boundary conditions is taken 
from the initial data and held constant. The fully sec- 
ond order characteristic variables always contain terms 
like dof for the evolved variables, and may be solved for 
evolution equations for (a, /3 r , x, 7rr, It) at the boundary 
if one additionally adds the D constraint to the bound- 
ary conditions. We could use the same method to derive 
constraint preserving conditions, but for now we hnd that 
the maximally dissipative conditions are satisfactory. At 
typical resolutions we find that the incoming constraint 
violation is of the order 10~ 7 — 10 -8 (see Fig. (j4])). 

Hydrodynamics equations are solved with an High- 
Resolution-Shock-Capturing (HRSC) method based on 
cell-center discretization and on the Local-Lax- Friedrichs 
(LLF) central scheme for the fluxes. The method is 
described in [3(1 HJ and we refer the reader to these 
references for all of the details. The assessment of the 
LLF flux in case of neutron star evolutions in full GR 
is presented in [32j . Different reconstruction schemes 
were tested, overall: the linear TVD based on Van- 
Leer Monotonized Centered limiter (MC2), the Convex- 
Essentially-Non-Oscillatory 3rd-order method (CEN03) 
and the Piecewise Parabolic Method (PPM). No relevant 
differences were found and the data presented were ob- 
tained with the CEN03 reconstruction. The metric fields 
are reconstructed with either simple averages or upwinded 
Lagrangian quartic reconstruction. 

The vacuum treatment is done with the use of an ar- 
tificial atmosphere as described in (29|. When, during 
the recovery of primitive variables, a point finish be- 
low a certain threshold pthr = /thr max (p(t = 0)) all 
the variables are set to the atmosphere level pi cv — 
j 'lev max (p(t = 0)) < pthr- Typical values are listed in 
Tab.|U 

In the collapse simulations matter fields were set to the 
atmosphere value once the horizon is formed and if a < 
^hydro-ex = 0.02. This is necessary to avoid numerical 
problems with the first points close to the origin. 

The code has been extensively tested for different prob- 
lems, especially the HRSC scheme. As expected, 4th 
order accuracy is reached in vacuum, while matter sim- 
ulations are accurate at 2nd order (except for situation 
with shocks). 

Stellar initial data are computed with a 2-domain 
pseudo-spectral method to solve the equations of Ref. [22| 
in isotropic coordinates, in spherical symmetry. A T = 2 
polytropic EoS: 

P{p) = (r - l)pe with e = K± - (B3) 

is used with K = 100 which is compatible with the 
ideal gas EoS used for the evolution. Chebyshev-Gauss- 
Lobatto grids with 64 collocation points per domain are 
sufficient for having the equations satisfied to machine 
accuracy. 
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